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A NON-DISSIPATIVE STAGGERED FOURTH-ORDER ACCURATE EXPLICIT FINITE 
DIFFERENCE SCHEME FOR THE TIME-DOMAIN MAXWELL’S EQUATIONS* 

AMIR YEFETt AND PETER G. PETROPOULOS 


Abstract. We consider a divergence- free non-dissipative fourth-order explicit staggered finite difference 
scheme for the hyperbolic Maxwell’s equations. Special one-sided difference operators are derived in order to 
implement the scheme near metal boundaries and dielectric interfaces. Numerical results show the scheme 
is long-time stable, and is fourth-order convergent over complex domains that include dielectric interfaces 
and perfectly conducting surfaces. We also examine the scheme’s behavior near metal surfaces that are not 
alligned with the grid axes, and compare its accuracy to that obtained by the Yee scheme. 

Key words. Maxwell’s equations, staggered schemes, finite differences, FD-TD scheme, explicit fourth- 
order schemes. 

Subject classification. Applied Mathematics 

1. Introduction. Recent engineering advances have resulted in ultra- wideband electromagnetic sources 
that find application in pulsed radar devices, ground-penetrating imaging systems, non-destructive evaluation 
of concrete structures, electronic on-chip interconnects, and novel communication systems. The need to 
simulate such problems requires fast and accurate solvers of the time-domain Maxwell’s equations in complex 
open/closed domains filled with heterogeneous dielectrics in which metals may be embedded. A mini-review 
of how the Computational Electromagnetics (CEM) state of the art impacts such technologies can be found 
in [13]. For about a decade, Yee’s Finite-Difference Time-Domain (FD-TD) algorithm [1] has provided the 
best [15] second-order accurate non-dissipative direct solution of the time-domain Maxwell’s equations on a 
staggered grid. The numerical error is controlled solely by the mesh size, and the algorithm is particularly 
easy to implement in the presence of heterogeneous dielectrics and metal boundaries. As with all finite 
difference schemes, the algorithm is inherently dispersive and anisotropic [4] and, for large-scale problems or 
for problems requiring long-time integration of Maxwell’s equations, errors from dispersion and anisotropy 
are significant unless the spatial discretization is extremelly small [10]. This leads to prohibitive memory 
requirements and high computational cost when addressing realistic problems. For some time now r , workers 
in CEM have realized the promise of high-order finite difference schemes. The question of staggered vs. 
collocated high-order schemes has been studied in [5] where it w r as shown that staggering is more efficient. 
At the same time a staggered high-order method can be constructed by altering a code that implements 
the staggered second-order accurate FD-TD algorithm. How'ever, the extended spatial stencil of high-order 
methods has inhibited their wide acceptance as it does not allow for easy application of boundary conditions 
(absorbing, metal) and modeling of dielectric interfaces. The issue of implementing an absorbing layer in 
a promising staggered scheme that is fourth-order accurate in space and second-order accurate in time [12] 
has been addressed in [14]. It remained though that this particular high-order scheme (desirable due to its 
similarities to the Yee scheme) was only second-order convergent and slightly more accurate than FD-TD 
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when used to solve problems with metal boundaries and/or dielectric interfaces. In this paper we revisit the 
explicit non-dissipative staggered scheme presented in [12] in order to address the remaining objections to 
its use. We demonstrate that it is divergence- free, and propose a series of numerical boundary conditions, 
involving one-sided differentiation and extrapolation, to implement metal boundaries and dielectric interfaces. 
The treatment of dielectric interfaces herein is different from that in [11] where a simple pointwise specification 
of dielectric properties was used; we show in our Numerical Results section that such an approach degrades 
the global convergence rate of the scheme to second order. Also, the treatment of metal boundaries herein 
is different than that used in [10] where the method of images was applicable due to the infinite extend of 
those boundaries in the numerical test performed there. As shown in our Numerical Results section, these 
new numerical boundary conditions give almost the same results as the compact-implicit Ty(2,4) scheme 
[2], and (more importantly) preserve the fourth-order accuracy of the scheme. Together with the results on 
absorbing layers [14], the present paper should open the way towards a general acceptance of this scheme, 
which is an extension of Yee’s algorithm to fourth-order accuracy. We now briefly outline the remainder of 
our paper. Section 2 describes the system of partial differential equations for which we will present the new 
boundary treatment for the fourth-order scheme considered herein. In Section 3 we present the details of 
the scheme and its numerical stability and accuracy properties; a derivation of the divergence-free property 
of the scheme is also given. Extensive numerical tests are given in Section 4; closed and open problems are 
considered and the actual convergence rate of a code based on the work herein is determined. In Section 5, 
we give a computational cost comparison between our approach and those of [1] and [2], Section 6 doses 
the paper with a short discussion, conclusions drawn from the numerical experiments, and a description of 
future work required to turn this scheme into an engineering tool. 

2. Preliminaries, The Maxwell equations in an isotropic non-dispersive medium are: 

dB 

V x E + — — 0 (Faraday's Law), 

( 2 . 1 ) 9t 

dD 

-Vx/fl = 0 (Ampere's Law), 

B = fiH, 

D = eE , 

coupled with Gauss’s law 

VB = 0 
V • D = p. 

To simplify the notation we shall consider the two dimensional case, with the only sources for the problem 
being incident waves. These waves will be scattered after they encounter an obstacle. Furthermore, in free 
space, c and // are constants. The extension of the system of equations to three space dimensions, and the 
inclusion of sources and variable coefficients is straightforward. In two dimensions, the system (2.1) now 
decouples into two independent sets of equations. We shall consider the Transverse Magnetic (TM) set of 
equations where the electric field is a scalar while the magnetic field is a two-dimensional vector. Letting 
r = ct = and Z = where e and fi are the permittivity and permeability coefficients in free space, 

and c is the speed of light, the TM equations are: 

0E Z dH y dH x 
dr 1 dx dy } 
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( 2 . 2 ) 


dHy _ 1 dE . 
dr Z dx 

3. The Scheme. The Yee scheme applied to (2.3) is 


e:x) = Kij + z ^ h :1Y 2 - Y 2 

Tjn+1/2 _ rjn-]/2 _ rn 

n x,iJ-l/2 ~ n x,i,j- 1/2 z Ay " 


H n+l { [% . = 

y,i— 1/2, j y, 1-1/2 


, _Li^La T? 71 
•1 + Z & X S * E *'*-W' 


where 


^ i J ^ *+ 1 /2 , j ' & i - 1 /2 , j 

(3.2) dyLjj = t^ij + 1/2 — Coj+1/2- 

U is evaluated at the appropriate time and space location. In order to improve the accuracy of this scheme we 

replace the difference operators (3.2) by the following fourth-order accurate stencil for the spatial derivatives, 

e£ 1- 
^■b*? d y- 

explicit(2,4): 


dU 

dy iJ+1/2 


(Uij-i - 27 Uij + 27^ J+ 1 - Uij+ 2 ) 


Hereafter we shall refer to (3.2) as the Yee scheme and to (3.3) as the explicit(2,4) scheme. To complete the 
fourth-order scheme at the first and last points of a bounded spatial domain we use fourth- and fifth-order 
accurate one sided appropximations to the derivative. We note that this is used only in order to globally 
approximate the derivative. No physical boundary conditions are included at this stage. These one-sided 
approximations are as follows: 


dy i,i/2 24Ay 


(— 22f7t,o + 17Ui t \ + 9U ia ~ 5^,3 + Co , 4 


1 = (- 23£/ M/2 + 210, 3/2 + 30.5/2 - 0,7/2) 


dU _ 1 

% i,N - 1 24 ^y 

dU_ _ 1 

dy j,jv— 1/2 ” 


(230,a'-i/ 2 ~ 210, /v— 3/2 — 30 ,v- 5/2 + 0,A r — 7/2) 


(220.1V — 170.1V-1 — 90.AT-2 + oUi'N-3 — Uj'N- 4 ) 


We next define 


Ah — — 
24 


-23 21 3 -1 

1 -27 27 -1 

0 1 -27 27 -1 


1 -27 27 -1 

1 -3 -21 23 
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(3.5) 


Ae =24 


-22 

17 

9 

-5 

1 


0 

1 

-27 

27 

-1 



0 

0 

1 

-27 

27 

-1 


0 

0 



1 

-27 

27 

-1 

0 


-1 

5 

-9 

-17 

22 


so the matrix form of the approximations to the derivative at the midpoint between grid points, and at the 
grid points is respectively: 



Ui/t 


f() 

d 

u 3/2 

_ 1 Ap 

U x 

dy 

• 

24Ay E 

Ujsi—i 


. OJV-1/2 


In 


0_ 

dy 


With these definitions, the matrix form of the discrete TM equations is: 


Ui 


C I/2 

U 2 

= 245„ Ah 

^3/2 





— 3/2 

Cat-i 


_ Un- 1/2 _ 


\EZu] n+l = [EZi,j] n + A „[HY i+a/2J ]" +1 /2 _ | [H X ij+1/2 ]”+V2 A i 


[HX iJ+l/2 ] n+l ' 2 = [HXij+w]*- 1 '* - —[EZi^A.^ 

[HY i+l/2 ,j}” +1 / 2 = [HY i+i/2J r-'/ 2 - A E [EZij] n . 

Note, the staggering in time and space results in a scheme that is second-order accurate in time and fourth- 
order accurate in space. — 

3.1. Divergence of Computed Fields. We now demonstrate that the explicit(2,4) scheme is divergence- 
free for TM waves, that is 

f r 0 * ' [ 

—div{iiH x ,nHy) = 0. 

Wo note that // may be a function of the spatial variables and introduce the relations: 

[{nHX) id+l/2 } n+ 'l* = [(nHX) iJ+l/2 ] n -^ - ^[EZ^r Afe 
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[(tiHY) i+l/2 j] n+l ' 2 = [(nm-) i+1/2J ] n - 1 ' 2 - — A e [EZ,j]". 

Multiplying the first equation above by 1 /Ax and A^, the second equation by 1 /Ay and Ag, arid adding 
we get 

(^A E [0,HX) ij+1/2 ] + ^[G<HY) i+1/2j ]A t E )"+ 1 = 


(^A E [(/,HX) iJ+1/2 ] + ^[(/^HY) i+1/2J ]A‘ ; ) n . 

Hence, if the field is numerically divergence-free initially it will remain so ever after. If the // is discontinous 
then the derivation of the property differs. In that case we segment the domain into subdomains, and on 
the boundaries of the dielectric we use a fourth-order extension of the approach developed for second-order 
schemes in [6] and [7]. The details of the fourth-order extension to handle dielectric interfaces is given in 
Section 4. 


3.2. Time Step. For each of the methods described above one must choose a time step for the numerical 
integration. This time step is based on two considerations: stability and accuracy. Since all the schemes 
have stability limits, this places an upper bound on the usable time step. 

The amplification factor is given by: 

& leapfrog = 1 ~h ~Z i )J ( 1 + — £ 2 ) 2 ~ 1 

r , , d/2 

z = A At. A is an eigenvalue of the Fourier transform of the spatial approximation. Let D = 

Then 

IX , _ 1 , [27sm(f ) - sin{f )] 2 [27sin(f ) - sin(f)] 2 1/2 ^7 

Inexplicit] 12 v ( Ax )2 + (Aj/) 2 ’ "3 

Since we wish to obtain higher order accuracy it is also necessary to limit the time-step by accuracy 
requirements. We do not want the accuracy of the scheme to be determined by the time integration. Hence, 
the temporal error should be equal to or less than the spatial error. For the Yee scheme one should choose 
the time step close to the stability limit. For the explicit(2,4) scheme the time step chosen depends on 
the accuracy desired. As the mesh is refined the spatial error decreases as of a fourth order scheme and 
so decreases faster than the temporal error. Thus, the time step should depend on (Ax) 2 . If the error 
requirements are too severe then this is inefficient and the leapfrog in time should be replaced by a fourth 
order Runge-Kutta method. However, for the experiments in this paper we shall use the same leapfrog 
method for both schemes. 


3.3. Spatial Accuracy. There are several ways of constructing fourth order methods. We can use 
either a staggered mesh or locate all the variables at the same mesh point. In addition we can either use a 
five point stencil in each direction to approximate the first derivative or else use a three point stencil with 
an implicit matrix inversion. Define the following operators: 

# DqU = VH/2p--l/2 
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« D >i — ^ U ^-3/2-^^3/2) + 27 {» j + i / 2-U j -1/2) 

• u i a — 04 h 

• iP?“li±} . ±L D .**h-' + ii(D 4 u)j = Si±U2ZSi=ui 

We define the truncation error as T = Du — = t / i 4 . By a Taylor series expansion we get 

• r i = efo ~ -004C875 

• r 2 = 5760 ‘ ^ 00295 

Kreiss and Oliger [3] give a simple analysis, for semi-discrete approximations, to calculate the number 
of points per wavelength needed to resolve a wave with speed a and given accuracy f. Let uj be the wave 
number (u(x,t) — e lu;x ~ af ). We consider the solution over q periods so that t = “2. They find that the 
number of points per wavelenth, A/, needed for accuracy e is given by 


M 


2n(2nr) l ^ p j 


Q\ 1/p 


where p is the order of the scheme. 

If we choose e = .01, one percent error, then we get 

• M 0 ~32.15<? 1/2 

• &Ii - 8.23r// 4 

• A/ 2 - 7.33 q ]/4 

As seen from the formulas of Kreiss and Oliger the benefits of a fourth order method, compared with a second 
order method, improve as one demands higher accuracy (e.g. 0.1 %) and with longer times of integration. 
We conclude that either staggering or a compact implicit method gives substantial improvement over the 
simplest fourth order accurate method. Combining staggering with an implicit method gives a little more 
improvement (see also [5]). We stress that staggering also helps in the imposition of the boundary conditions. 
In addition using the Yee placement of the variables simplifies the conversion of an existing code to fourth 
order accuracy. 

If one assumes a uniform grid-spacing, i.e, Ax — Ay — h, and defines the number of points per wavelength 
to be N = then the numerical wave speed c* can be written as: 


r * = — 

C Y ee _ 

7 T 

t-cxphcit 24 tt 

+ 


. 2 , 7T COs((?) . 2, 

sm 2 (— — T2) + sm 2 ( 


N 


7rsin(0) 1 

N ' 


. ^cosfflK . , 37rcos(0Y 
27 sin ( — ~ sm( 


N 


. ,7rsin(6») . 3ttcos(i9) 

27 S1I1 ( 5T7 ) - Sm ( AT ) 


N 

l2 


N 


N 


B/ 2 


Shown in Fig. 3.1 are the polar diagrams of the numerical phase speed for A r = 1,2, 4, 8 and 16 for the two 
schemes. A comparison of the numerical phase speed for N = 20 is given in Fig. 3.2. The numerical phase 
speed in the two schemes experiences a phase lag. The lag decreases §s N increases. The lag decreases for 
the explicit scheme. In other words, for a given grid-spacing, the error (1 — c*) for high frequency modes is 
greater than that for low frequency modes. For a fixed N, the error using Yee’s scheme is greater than The 
explicit. The two figures also demonstrate the anisotropy inherent in the discretizations. 

One observes that the error is the greatest along the axes ( 6 = 0, tt/2, 7 r and 37r/2) and the least along the 
diagonals (9 — 7t/4,37t/4, 57t/ 4 and 77r/4). An important quantity to measure is the isotropy error defined 
as the difference of the maximum and the minimum values of the numerical phase speed. For N— 20, as 
shown in Fig. 3.2 the isotropy errors are 0.2% for the Yee scheme and 0.0034% in for the explicit scheme. 
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Expanding c* one obtains the error: 

1 - C Yee = (g + 24 COS ^^ 

The above equation shows that the leading dispersive errors in Yee’s scheme is inversely proportional to N 2 
and in the explicit scheme is inversely proportional to N 4 , This, of course, is just a reflection that Yee’s 
scheme is second order accurate in space while the explicit scheme is fourth order accurate in space. It also 
shows that the leading dispersive error is a function of the wave propagation direction 0 on each grid. 



Fig. 3.1. Polar diagram of the numerical phase speed. Yee: , explicit (2 ,^):— • 



Fig. 3.2. Comparison of numerical phase speed } 0 < 9 < 2n. 

3.4. Numerical Errors of Spatial and Time Discretization. We have discussed spatial and time 
errors individually in the previous sections. In this section we will investigate the numerical errors from the 
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combined discretization. These errors can be determined from the eigenvalues of the amplification matrix. 
In Fig. 3.3 the curves labeled ’exact 5 are the normalized phase shift using the exact time integration. The 
difference between these curves and the value 1 is the spatial phase error, which can have either a phase lag 
or a phase lead. The other curves are the normalized phase shift using the staggered leapfrog (with Yee and 
the explicit(2,4)) with CFL — For the Yee scheme the difference between each curve and the exact curve 
is an additional phase error due to time discretization. Leapfrog has a larger error and a phase lead and so 
moves the curve away from 1. Except for a very small contribution due to the nonlinear relation between a 
and A, the time discretization does not reduce the isotropy error introduced by the spatial discretization. 



Fig. 3.3. Comparison of phase shifts for complete discretization. 


4. Computational Results. In this section we compare three different schemes: the standard Yee 
scheme, the explicit,(2.4) scheme as extended herein, and the Ty(2,4) [2] scheme. All schemes are advanced 
in time by the leapfrog method. For all computations we choose Z = 1. For the numerical examples that are 
posed in an open domain we use a PML method in the far field. In all cases the error is measured against 
the exact E z in the L > norm. All numerical examples herein test the accuracy and stability of the one-sided 
differencing and extrapolation operators introduced in Section 2 and in the present Section. 

We first consider a test case with the following initial and boundary conditions: 

E z = sin(37rx) sin(47 xy) 

H y = (3/5) cos(37rx) sin(47n/) sin(^^) 

7T /\ t 

H x = -(4/5)sin(3^x)cos(47ry)sin(— — ) 

The exact solution in this case is: 

E z = sin(37rx) sin(47n/) cos(o7rt) 

For the three schemes we choose uniform grid spacing. For the Yee, the explicit(2,4) and the Ty(2,4) schemes 
we take: h — Ax = Ay = 1/20, 1/40, 1/80. For the Yee scheme we set At = 2h/3, while for the explicit (2,4) 
and Ty(2,4) schemes we set At = h 2 , Figures 4.1a-4.1c show the actual logarithmic errors as a function of 
time. Figure 4. Id shows the convergence rate of the L 2 spatial error in the maximum norm over the time 
interval [0, 10]. The slope of the line for the Yee scheme converges to 2, and the slope of the lines for both 
the explicit(2,4) and Tv(2,4) schemes converge to 4. This can also be seen in Table 4.1. 
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TIME TIME 


Fig. 4.1a. Iog io (||err 0 r||£ 2 ) For the Yee scheme. Fig. 4.1b. Iog lo (||err 0 r||£ 2 ) for the explicit(2 f 4) 

scheme . 



Fig. 4.1c. log lo (||e7T0r||L 2 ) for the Ty(2,4) scheme. Fig. 4. Id. Logio(\\error\\i 2 ) as a function of 

Logio(h) For the Yee, the explicit(2,4) and theTy(2,4) 
schemes. 


Table 4.1 

The maximal errors in L 2 norm for each mesh size 


scheme 

h 

At 

Max(\\error\\L 2 ) 

rate 




0 < t < 10 


explicit( 2, 4) 

1 

20 

1 

400 

0.014 


explidt( 2, 4) 

1 

40 

1 

1600 

1.43 x lO” 4 

6.60 

explicit( 2, 4) 

1 

80 

1 

3200 

4.76 x 10“ 6 

4.9152 

Ty( 2,4) 

1 

20 

1 

400 

0.0224 


Ty( 2,4) 

1 

40 

1 

1440 

7.07 x 10~ 5 

8.306 

Ty( 2,4) 

1 

80 

1 

3440 

2.12 x 10“ 6 

5.0589 

Yee 

1 

20 

1 

35 

0.189 


Yee 

1 

40 

1 

70 

0.0475 

1.99 

Yee 

1 

80 

1 

140 

0.0118 

2.003 


We next consider the simplest mode of propagation in a rectangular cross section wave guide. We take 
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the walls to be perfect conductors. We take the following boundary and initial conditions: 


E z (x,y>0) = sin(37rj*) sin(47ry) 

rr , At 3 . 57rAt . /4 , 

H u {x,y,—) = --sin(3jrx — )sin(4 iry) 

rr / At, 4 , 5 ttAT . , 

H x {x,y, — ) = - - cos(3;rx — ) sm(4ny) 

£~(0,i/, f) = — sin(57rt) sin(47 ry) 

E z (l,y, t) = sin(37r — 5 t rt) sin(47ry) 

E z (x > 0, t) = 0 

E z (x, M) = 0 


The solution is then: 


E z (x , y , t) = sin(37rj: — 57rt) sin(47 ry). 


For the three schemes we choose uniform grid spacing. For the Yee, the explicit(2,4)> and the Ty(2,4) schemes 
we again take: h = Ax = A y = 1/20, 1/40, 1/80. The CFL numbers are chosen as before. In the figure 4. 2d 
we draw the error as a function of the mesh size. For this test problem too the slope for the Yee scheme 
converges to 2 while that for the explicit(2,4) converges to 4. In table 4.2 we present the errors in L 2 norm 
for the Yee and the explicit (2,4) scheme. In both cases the errors are almost linear in time. However as the 
mesh is refined the Yee scheme yields second order accuracy while the explicit (2,4) yields between fourth 
and fifth order accuracy. 




Fig. 4.2a. log l 0 (||error||x. 2 ) F° r Fee scheme. Fig. 4.2b. log 10 (|[error||^ 2 ) For the explicit(2,4) 

scheme. 
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Fig. 4.2c. log l0 (ll error llL 2 ) For tk e Ty(2,4) 


Fig. 4. 2d. Logio(\\e.rror \\^ 2 ) as a function of 


scheme . 


Lo<7io(/ 0 For Zfte Yee, the explicit(2,4) and the.Ty(2,4) 


schemes . 


Table 4.2 

77ie maximal errors in L -2 norm. 


scheme 

h 

At 

M ax{\\error\\i 2 ) 

rate 




0 < t < 10 


explicit (2, 4) 

i 

20 

I 

400 

0.014 


explicit^ 2,4) 

1 

40 

1 

1600 

1.9316 x 10“ 4 

6.2 

explicit (2, 4) 

1 

80 

1 

3200 

6.48 x 10 -6 

4.896 

Ty( 2,4) 

1 

20 

1 

400 

0.0242 


Ty( 2,4) 

1 

40 

1 

1440 

7.9304 x 10“ 5 

8.15 

Ty(2, 4) 

I 

80 

1 

3440 

2.329 x 10“ 6 

5.089 

Vee 

1 

20 

30 

0.1889 


Fee 

1 

40 

1 

60 

0.0476 

1.9885 

yee 

1 

80 

120 

0.0119 

2.0032 


We next consider the treatment, of a domain which contains air and a lossless dielectric with a relative 
permittivity of ti as shown in Figure 4.3. 



Fig. 4.3. The computational domain 
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Since E z is continuous while its second derivatives are discontinuous we use the following fourth-order 
explicit interpolation to implement the discontinuous dielectric properties in the explicit (2,4) scheme. Define 


(4.1) 


Then 


4. . — _ 
* tni ~ 16 


5 

-1 

0 

0 

0 


15 

9 

-1 


-5 

9 

9 


1 


-1 

1 


0 

-1 . 

9 9 

-5 15 


<4 


e l/2 

C‘2 

— 4 

f 3/2 


— ^int 


Cp-2 



. <*-1 - 


. e p- 1/2 . 


0 

0 

0 

-1 

5 


In two dimensions this is replaced by: 

[ € ij] ~ 2 + [ e ij+l/2]A t int) * 

For the Yee scheme we take ^interface = An exact solution in this case is: 


E~ = 


2cos(^-X) cos(ut) sin(KyY) | JV| <\ 0 < Y < 1 

exp(^S)exp(-2a^2|A'|) cos(u rf) am(K p Y) |A'| > \ 0 < Y < 1 


Hy= < 


-\/t 2 — fi sin(^ L A') sin(u;t) sin(A' J/ V') 




| A' | < | 0 < Y < 1 

exp( 2L )p) exp(— X) sin(^.’/) sin(AyY) A" > ^ 0 < 1* < 1 


' 1 1 ex p(2Lyd) exp( sin(cjf) sin(Jvp5 r ) 


A: <-i 0 < r < 1 


h t = 


—\/(i + 3e 2 cos(^pA') sin(u;f) cos(A' v l') |A'| < ^ 0 < Y < 1 

ex p(S#) exp(-^j.Y|) sin(wt) cos(K y Y) |A'J > \ 0 < Y < 1 


where K u = . / t i +3t a and u> = ~ ,* n ■ ■ Here we choose f| = 1 and e 2 = 2,4. We use the same mesh sizes 

y 3 V f2— *1 3Vf2— eT 7 

as before. In Figures 4.4 and 4.5 we draw the errors as a function of the mesh size for various 62. 




Fig. 4.4. The maximal errors in L 2 norm as a Fig. 4.5. The maximal errors in L 2 norm as a 

function of the mesh size with €2 — 2. function of the mesh size with C 2 = 4. 
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Although we get only second order accuracy for the explicit(2,4) scheme we still get much better results 


Tabi.f. 4.3 

The maximal errors in L 2 norm with C 2 = 2. 


scheme 

h 

At 

Max(\\error\\L 2 ) 

rate 




0 < t < 10 


explicit ( 2, 4) 

1 

20 

1 

400 

0.0019 


e:rpZzczf(2, 4) 

1 

40 

1 

1600 

5.7585 x lO" 4 

1.715 

explicit (2, 4) 

1 

80 

1 

3200 

1.4909 x 10“ 4 

1.94 

7y(2,4) 

1 

20 

1 

400 

0.00196 


ry(2,4) 

1 

40 

1 

1600 

5.7721 x lO" 4 

1.763 

Ty(2,4) 

1 

80 

1 

6400 

1.4955 x 10~ 4 

1.948 

Yee 

1 

20 

1 

30 

0.0363 


Yee 

1 

40 

1 

60 

0.0089 

2.028 

Yee 

1 

80 

1 

120 

0.00222 

2.003 


Table 4.4 

The maximal errors in L 2 norm with f 2 1=1 4. 


scheme 

/i 

At 

Max(\\error\\L 2 ) 

rate 




0 < t < 10 


explicit( 2, 4) 

1 

20 

1 

400 

0.0014 


explicit (2, 4) 

1 

40 

1 

1600 

3.765 x 10~ 4 

1.894 

explicit ( 2, 4) 

1 

80 

1 

3200 

9.7748 x 10“ 5 

1.945 

Ty( 2,4) 

1 

20 

1 

400 

0.00139 


Ty( 2,4) 

1 

40 

1 

1600 

3.756 x 10~ 4 

1.887 

r»(2,4) 

1 

80 

1 

6400 

9.7579 x 10” 5 

1.944 

Yee 

l 

20 

1 

30 

0.0095 


Yee 

1 

40 

1 

60 

0.00237 

2.003 

Yee 

1 

80 

1 

120 

5.9442 x 10“ 4 

1.9953 


than the ones we get employing the Yee scheme. However, we are using a fourth-order scheme and the loss 
of two orders of convergence in the presence of heterogeneous dielectrics is undesirable. 

An innovative approach to handle heterogeneous dielectrics for the Yee scheme is presented in [6]-[7]. 
We have extended this approach to include heterogeneous piecewise constant dielectric properties in the 
explicit (2,4) scheme. As we see below, the fourth-order convergence is recovered globally. We divide the 
computational domain into three subdomains. Two contain air and the third one contains the lossless 
dielectric. On the interfaces both the electric and magnetic fields are approximated as follows. Suppose the 
interfaces are located at i = Ii and i — I 2 , and e = 62 for I\ < i < I 2 while e = e\ for i > I 2 and i < I\ . We 
approximate H y at i = /1 and i — 1 2 by using the following fifth-order extrapolation: 

rjn+l/2 315 rrn+1/2 _ 105 rrn+1/2 , 189 rrn+1/2 _ *5 rrn + 1/2 , 35 rrn+1/2 

viij ~ 128 yjl_1/2j 32 ^ 1 - 3 / 2 j “ r 54 v/j- 5 / 2 , j 32 yi 1 - 7 / 2 , > ^ 128 y/ i- 9 / 2 ^' 

rjn+l/2 315 +] / 2 _ 12£rrn + l/2 , ^9 l+1/2 _ 45 „ n+ i/2 3 ^ rrn+1/2 

128 V,2 + 1/a,i 32 y/ 2+3/2,i ^ 04 F/ a +5/2J 32 V /2 + 7/2.i 128 y/ 2+9/2.i 
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Hyi l" 1 / 2 can be extarpolated by using the points to the left of I\ , or by using the points to the right of 
To because H y is a continuous function. Once H y is approximated on the interface we approximate the 
x-derivative of H y using H yil J and H yi2 3 as follows: 


_®_zjn + 1/2 ^ 30‘2 ——H ——H +—H 

dx viij 105 ' 3 8 24 y/l+3/2,i 4 q /7 ^o+s/ 2.> “t" 45^0+7/2., 

d rrn+l/2 _352 ^rr , 21 tt 5_ir 

dx yi2 j 105 y/2 ' J ^ 8 y/ 2- l /2 24 y#a “ s/2,i 40 ^ /2_5/2,j 4C y/2 * 7/2J 


Since we discretize the time we want to lose as little information as possible. Therefore, we approximate 
/Cty 2 where moves more slowly, i.e. where e — c* 2 - Once ^ Z/^ 1 / 2 and are calculated we 

can evaluate EZf+j and J EZfJ'j the following way: 


77771 + 1 


= ^.J- 


At 


A 


/l.i-3/2 


- 2777 


■P/pj-l/2 


+ 27 H, 


X 0. J+l/2 


■H 


a ‘/i-i+3/ 


.) 



f 352 rrrH-1/2 _ 33 itti+I/ 2 , 33 un+1/2 

\105 y/lJ 8 y o+i/2,i 24 y o+ 3 /2,i 


^lirn+l/2 , A/7"+l/2 | 

40 y o+5/2j ^ 45 y/ L +7/2j y 


Since ^ and are approximated where c = €2 we set t — to at i = and i — I>. ^Hy^J 2 

is approximated by using points, which are on the right-hand side of I \ . This is done because the velocity 
of the waves is y/r which is smaller than the velocity on the left-hand side of I \ . 

In Figures 4.6a-4.6c we draw the error of the three schemes for various mesh sizes with 62 = 2 and in 
Fig 4.7a-4.7c we draw the error of the three schemes for various mesh sizes with €2 ~ 4. In Fig 4.6d ami 
4.7d we draw the errors as a function of the mesh size. The slope of the Yee scheme is 2 and the slope 
of the Ty(2,4) scheme and the explicit(2,4) converges to 4 as can also be seen in table 4.5 and table 4.6. 



Fig. 4.6a. log 10 (||error||£, 2 ) For the Yet scheme. 


Fig. 4.6b. Iog 10 (||errt>r||L 2 ) f or the explicit(2,4) 

scheme. 
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Fig. 4.6c. log 10 (||error||£, 2 ) for the Ty(2,4) scheme. Fig. 1.6d. Log\o{\\error\\i 2 ) as a function of 

Logio(h) For the Yec, the explicil(2,4) and theTy(2,4) 
schemes . 


Table 4.5 

The maximal errors in hi norm with ei — 2. 


scheme 

h 

At 

Max{\\error\\i 2 ) 

rate 




0 < t < 10 


explicit (2, 4) 

l 

20 

1 

400 

3.1829 x 10 -4 


explicit (2, 4) 

i 

40 

1 

1600 

4.9839 x 10“ 6 

5.996 

explicit( 2, 4) 

1 

80 

1 

3200 

2.6518 x 10“ 7 

4.23 

Ty{ 2,4) 

1 

20 

1 

400 

1.978 x 10“ 4 


Ty( 2,4) 

1 

40 

1 

1600 

2.3593 x 10“ 6 

6.389 

Ty(2,4) 

1 

80 

1 

6400 

4.4066 x lO" 7 

2.420 

Yee 

1 

20 

1 

30 

0.0363 


Yee 

1 

40 

1 

60 

0.0089 

2.028 

Yec 

80 

1 

120 

0.00222 

2.003 




Fig. 4.7a. log 10 (||error||/ j2 ) For the Yee scheme. 


Fig. 4.7b. log 10 (||error[|£ 2 ) for the explicit(2,4) 
scheme. 
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Fig. 4.7c. Iog 10 (||error||^ 2 ) for the Ty(2 1 4) scheme. Fig. 4.7d. Logio(\\error\\i 2 ) as a function of 

Logio(h) For the Yee, the explicit(2,4) and theTy(2,4) 
schemes. 


Table 4.6 

The maximal errors in L 2 norm with c 2 = 4. 


scheme 

h 

At 

Max(\\error\\i 2 ) 

rate 




0 < t < 10 


explicit (2, 4) 

1 

20 

1 

400 

6.9239 x Hr 5 


explicit (2, 4) 

1 

40 

1 

1600 

3.5486 x 10~ 6 

4.286 

explicit (2, 4) 

1 

80 

1 

3200 

2.0112 x 10“ 7 

4.141 

Ty{ 2,4) 

1 

20 

1 

400 

2.7043 x 10“ 5 


Ty{ 2,4) 

1 

40 

1600 

1.4233 x 10“ 6 

4.249 

T 2/(2, 4) 

1 

80 

6400 

1.1040 x 10“ 7 

3.688 

Yee 

! 1 
20 

1 

30 

0.0095 


Yee 

1 1 
40 

1 

60 

0.00237 

2.003 

Yee 

1 

80 

1 

120 

5.9442 x 10“ 4 

1.9953 


Next, we look at a coated perfect conductor with a dielectric layer. The coating thickness is | with a 
relative permittivity of e 2 as shown in Fig 4.8 . 



Fig. 4.8. The computational domain 
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1. An exact solution in this case can be: 


We take c 2 = 2 and f\ 
E z 


sin(aiA”) sin(u;t) sin(6Y) 0 < X < ~ 0 < Y < 1 

cos(a >X) sin(a;£) sin(61 r ) \ < X < | 0 < V < 1 


where 


We choose: 



- 2J- cos(a] .V) cos (ujt) sin(fcY) 
2* cos(a 2 -V) cos(wf) sin(6Y) 

~ sin (a i A") cos(u;f) cos(6Y) 
j sin(a -2 A r ) cos(uj£) cos(6Y) 


0 < X < ~ 0 < r < 1 
| < a < f 0<r < 1 

0 < A < I 0 < Y < 1 
| < A r < f 0 < Y < 1 


aj + b 2 = e 2 u ) 2 

a? 2 + b 2 ~ 6ju ; 2 

x = ^ : sin(^) = cos(^) 

x=f: cos(^) = 0 


Ci — 1 
6 2 =2 
a\ — 37r 
02 = 27T 
b — Ti 

u — Vbn 


On the interface between the air and the dielectric we use the same technique we used in the previous ex- 
ample. In Figures 4.9a-4.9c we draw the errors for various mesh sizes. In Figure 4.9d we draw the errors as 
a function of the mesh size. For the Yee scheme we get second order accurcy and for the Ty(2,4) and the 
explicit (2,4) the accurcy converges to 4, which can also be seen in Table 4.7. 



Fig. 4.9a. Iog 10 (||error||£ 2 ) For the Yee scheme. 
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Fig. 4.9c. Iog 10 (||error]|z, 2 ) for the Ty(2,4) scheme. Fig. 4.9d. Log\o(\\error\\i 2 ) as a function of 

Logxo(h) For the Yee, the cxplicit(2,4) and theTy(2 f 4) 
schemes. 


Table 4.7 

The maximal errors in L<i norm . 


scheme 

h 

At 

A/aa:(||erro7’||i 2 ) 

rate 




0 < t < 10 


explicit (2, 4) 

1 

20 

1 

400 

0.00398 


explicit (2, 4) 

1 

40 

1600 

2.4868 x 10“ 4 

4.001 

explicit^ 2, 4) 

1 

80 

3200 

1.0889 x 10“ 5 

4.513 

Ty{ 2,4) 

1 

20 

1 

400 

0.0063 


Ty( 2,4) 

1 

40 

1 

1600 

1.609 x 10“ 4 

5.308 

Ty( 2,4) 

1 

80 

1 

6400 

1.7820 x 10“ 6 

6.497 

Yee 

1 

20 

1 

30 

0.1498 


Yee. 

1 

40 

60 

0.037 

2.004 

Yee 

1 

80 

1 

120 

0.0093 

2.0026 


To our knowledge, tills is the first fourth-order scheme that preserves its convergence rate when discon- 
tinuities in the coefficients are present. 

Next we consider a monochromatic isotropic point source of wavelength 0.25, that is switched on at 
t — 0 and radiates in frce-space. The domain is 0 < x,y < 1. For Yec’s scheme we choose h = -h, A f = ■—, 
for the explicit(2,4) and t.heTy(2,4) scheme h= A t = h 2 . The point source is modeled by adding a term 

representing a current I z (t) = 0.01 sin(8irt)e(t) at (a:,y) = (p |) where e(t) denotes the Heaviside unit-step 
function. For the pulse under consideration, the radiated field is the solution of 

(4-2) &*E Z + d]E z - d\E z = Zd t I z {t)S{r - r s ) 

r — r s = (x — — I). The solution consists of rotat tonally symmetric outgoing waves. The exact solution 

is 


E z {r,t) = 


Z_ 

2i r 


f 

Jo 


dthjt - \/Kr ~r s | 2 + £ 2 )) 




+ e 




In Fig. 4.10 we plot the errors, in the Li norm, for the various approximate solutions. 


18 




x 10' 



TIME 


Fig. 4.10. The errors in L 2 norm between the numerical solution and the exact solution of 4-2. 



Fig. 4.11. The contours of the Yee scheme and 
the exact solution at t = 8. The Yee scheme is drawn 
by and the exact solution is drawn by a solid line 


Fig. 4.12. The contours of the explicit(2 r 4) 
scheme and the exact solution at t = 8. The ex- 
plicit(2,4) scheme is drawn by and the exact so- 

lution is drawn by a solid line. 


Next we consider a monochromatic isotropic point source of wavelength 0.25, that is switched on at t = 0 
in the presence of an infinite perfect conductor( fig 4.13). The domain is —00 < x,y < ^ —00 < y < 00 . For 
Yee’s scheme we choose h = for the explicit(2,4) and theTy(2,4) schemes h = At = h 2 . The 

point source is modeled by adding a term representing a current I z (t) = 0.01 sin(87r t)e(t) at {x t y) = j) 
where e(t) denotes the Heaviside unit-step function. The equations in this case arc: 


0E* z( dHy 

Ot [ dx 
dH x _ 1 dE z 

dt Z dy 
dHy _ 1 dE z 
dt Z dx 


dH x 

dy 


) - ZI z (t)6(r - r 8 ) 


With the boundary conditions: 

(4.3) E z (l/2, y,t) = 0 
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FlG. 4.13. The computational domain. 


The exact solution can be constructed by using the exact solution for the previous case and by using the 
image method. 


In Figures 4.14, 4.15 and 4.16 we draw the errors in L 2 norm and the contours of the exact and numerical 
solutions. 



Fig. 4.14. The errors in L 2 norm between the numerical solutions and the exact solution of J t .3. 
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Fig. 4.15. The contours of the Yee scheme and Fig. 4.16. The contours of the explicit(2,4) 

the exact solution at t = 5. The Yee scheme is drawn scheme and the exact solution at t — 5. The ex- 

by and the exact solution is drawn by a solid line plicit(2 r 4) scheme is drawn by and the exact so- 

lution is drawn by a solid line . 

Next we consider a monochromatic pointsource in the presence of an inclined perfect conductor Fig 4.17a 
. We test these three schemes by using the staircasing method. We meshure the error in L\ norm along the 
dashed line which can be seen in Fig 4.17b . As can be seen from Fig 4.17c the Ty(2,4) scheme is unstable 
whereas the Yee scheme and the explicit.(2,4) schemes are stable. In [9] the Ty(2,4) scheme was tested as 
w r ell but there the Ty(2,4) scheme w r as stable. 



Fig. 4.17c. The [ |error [| X, i as a function of time for the Yee, the explicit(2,4) and theTy(2,4) schemes . 
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Next we test these three schemes in [0, 1/2] x [0, 1/4] x [01/2] domain. An exact solution in this case 
can be [8]: 

H x 
Hy 

H z 

E x 

Ey 

E z 

Where 

J 1 = .4 2 + B 2 + C 2 

0=A+B+C 

We choose: 

,4 = 7 r 
B = -2t r 
C — 7T 

id = s/Gtt 

In table 4.8, we can see that for the explicit(2,4) and the Ty(2,4) schemes we have used At — h 2 and for 
the Yee scheme we have used At = — . The explicit (2,4) as well as the Ty(2,4) schemes behave better than 
expected and gives almost fifth order of accuracy. 

For all this cases we have measured the error between the approximated electric field in the z direction 
and the exact electric field in the z direction in norm. 


Table 4.8 

Comparison of the errors in L2 norm, 


scheme 

/1 

At 

Max(\\error\\t 2 ) 

rate 




0 < * < 10 


explicit^ 2, 4) 

1 

20 

400 

5.375 X 10- 4 


explicit ( 2, 4) 

1 

40 

1600 

2.184 x 1(T 5 

4.621 

explicit( 2, 4) 

1 

80 

3200 

9.071 x 10~ 7 

4.590 

Ty{ 2,4) 

20 

1 

400 

3.621 x 10“ 4 


Ty( 2,4) 

1 

40 

1 

1600 

1.144 x 10“ 5 

4.983 

Ty(2,4) 

1 

80 

1 

6400 

3.5621 x 10“ 7 

5.005 

Fee 

1 

! 20 

1 

35 

0.0027 


Fee 

1 

40 

1 

70 

7.3 x 10“ 4 

1.9028 

Fee 

1 

80 

1 

140 

1.8252 x 10“ 4 

2.0015 


— sin(urt) sin(Ax + By -f Cz) 

— sir \(ut) sin (At + By + Cz) 

= sin(cjf) sin (Ax -F By + Cz) 

C — B 

— cos(u;f ) cos (Ax + By 4- Cz) 

U) 

A - C 

= cos(u;<) cos(At 4- By + Cz) 

LO 

B - A 

— cos(u;f) cos (At + By + Cz) 
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Fig. 4.18a. log^d (error ||j, 2 ) for the Vee scheme. 


Fig. 4.18b. logio(Il eT ’ ror lk 2 ) f or ^ e ex ~ 
plicit(2,4) scheme. 




Fig. 4.180. log 10 (|[error||^ 2 ) for the Ty(2,4) Fig. 4.18d. Log\^{\\error\\i 2 ) as a function 

scheme. of Log\o(h) For the Yee, the explicit(2,4) and the 

Ty(2,4) schemes. 

5. Computational Cost Comparisons. In order to compare the efficiency of the explicit (2, 4), the 
Ty(2,4) and Yee scheme we examine the following boundary conditions: 


E z = sin(37ra:) sin(47ry) 

H y = (3/5) cos(37nr) sin(47n/) sin( 


SitAt 


2 y 

H x — — (4/5) sin(37rj*) cos(47ry) sin( — - — ) 


The exact solution in this case is: 


E z = sin(37nr) sin(47 ry) cos(57r t) 

For the Ty(2,4) sclieme and the explicit (2,4) scheme we use a uniform gridspacing with Ax = Ay — -~ w For 
the Yee scheme we also use uniform gridspacing with Ax = Ay = —q. We chose these mesh sizes in order 
to get the same error between the exact E z and the approximated E z in L 2 norm. The comparison is shown 
in table 5.1. The programs were written in fortran and run on a Digital Alpha workstation. 

The CPU time needed to achieve the same accuracy in Yee’s case is more than 11 times larger than 
required for the Ty(2,4) scheme and 91 times larger than required for the explicit(2,4) scheme. 
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Table 5.1 

CPU -time using various difference schemes. 


scheme 

h 

At 

Max(\ | error | \i 2 ) 
0 < t < 10 

CPU — time 

explicit (2, 4) 

Ty[ 2,4) 
Yee 

i 

30 

1 

30 

1 

240 

1 

900 

1 

900 

1 

360 

i.99 x nr 3 
i.25 x nr 3 
1.31 x 10“ 3 

0.9 sec 
5.7 sec 
91 sec 


6. Discussion and Conclusion. The results demonstrate that we can use a coarser mesh with the 
fourth order scheme and still get the same accuracy as with the Yee scheme. This is true even in the presence 
of a dielectric media. 

Although this scheme is not as good as the Ty(2,4) schcme[2], it is still easier to modify an existing code 
based on the Yee scheme and make it fourth order accurate, by using the explicit(2,4) scheme. This is true 
because in the Ty(2,4) scheme one has to inverse a matrix by using a LU decomposition. 
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